home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_5 / a5_3a.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.9 KB  |  121 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 5.3 (Non-linear Curve Fitting).
  9. % A detailed look at the exponential curve fit.
  10. % Section    5.2,    Curve Fitting, Page 280
  11. echo on; clc; format short; hold off; clear
  12.  
  13. % This program performs exponential curve fitting,
  14.  
  15. % by using the method of data linearization. Given a set
  16.  
  17. % of data points { (x , y ), (x , y ) ,..., (x , y ) },
  18. %                    1   1     2   2          n   n
  19.  
  20. % The abscissas and ordinates are stored in X and Y, respectively.
  21.  
  22. % X = [x , x  ,..., x ];  Y = [y , y  ,..., y ];
  23. %       1   2        n          1   2        n
  24.  
  25. pause % Press any key to continue.
  26.  
  27. clc;
  28. %    Place the abscissas for the points in  X.
  29.  
  30. %    Place the ordinates for the points in  Y.
  31.  
  32. X = [0    1    2    3    4];
  33.  
  34. Y = [1.5  2.5  3.5  5.0  7.5];
  35.  
  36. pause    % Press any key to graph data points.
  37.  
  38. clc; clg;
  39. axis([-0.1 4.1 -0.2 8.2]);
  40. plot(X,Y,'or');...
  41. hold on;...
  42. plot([-10000 10000],[0 0],'b',[0 0],[-10000 10000],'b');...
  43. xlabel('x');...
  44. ylabel('y');...
  45. title('The given x-y data points');...
  46. grid;...
  47. axis;...
  48. hold off;...
  49. shg; pause    % Press any key to continue.
  50.  
  51. points = [X;Y];
  52. format short;
  53. clc,echo off,diary output,disp('The given x-y points:'),...
  54. disp('      x         y'),disp(points'),diary off,echo on
  55. pause    % Press any key to "linearize" the data points.
  56.  
  57. clc; clg;
  58. XT  = X;
  59. YT  = log(Y);
  60. Ya = 1;
  61. p = polyfit(XT,YT,1);
  62. A = p(1);
  63. B = p(2);
  64. Xs = [min(X) max(X)];
  65. Ys = polyval(p,Xs);
  66. axis([-0.2 4.2 -0.1 2.1]);...
  67. plot(XT,YT,'or',Xs,Ys,'-g');...
  68. hold on;...
  69. plot([-10000 10000],[0 0],'b',[0 0],[-10000 10000],'b');...
  70. xlabel('X');...
  71. ylabel('Y');...
  72. Mx1 = 'Least squares line: Y = ';...
  73. Mx2 = [Mx1,num2str(A),' X'];...
  74. if B > 0,
  75.   Mx2 = [Mx2,' + ',num2str(B)];
  76. else
  77.   Mx2 = [Mx2,' - ',num2str(abs(B))];
  78. end;...
  79. title(Mx2);...
  80. grid;...
  81. axis;...
  82. hold off;...
  83. shg; pause    % Press any key to continue.
  84.  
  85. points = [XT;YT];
  86. format short;
  87.  
  88. Mx1 = 'The transformed points in the X-Y plane.';
  89. Mx3 = '      X         Y';
  90. clc,echo off,diary output,...
  91. disp(''),disp(Mx1),disp(Mx2),...
  92. disp(Mx3),disp(points'),diary off, echo on
  93. pause    % Press any key to fit y = C exp(Ax).
  94.  
  95. clc; clg;
  96. A = p(1);
  97. B = p(2);
  98. C = exp(B);
  99. h = (max(X)-min(X))/150;
  100. Xs = min(X):h:max(X);
  101. Ys = C*exp(A*Xs);
  102. axis([-0.1 4.1 -0.2 8.2]);...
  103. plot(X,Y,'or',Xs,Ys,'-g');...
  104. hold on;...
  105. plot([-10000 10000],[0 0],'b',[0 0],[-10000 10000],'b');...
  106. xlabel('x');...
  107. ylabel('y');...
  108. Mx1 = 'The fit f(x) = ';...
  109. Mx2 = [Mx1,num2str(C),' exp(',num2str(A),' x)'];...
  110. title(Mx2);...
  111. grid;...
  112. axis;...
  113. hold off;...
  114. shg; pause    % Press any key to continue.
  115.  
  116. points = [X;Y;C*exp(A*X);Y-C.*exp(A*X)]';
  117. Mx3='    x(k)      y(k)      f(x(k))   error';
  118. clc,echo off,diary output,disp(''),disp(Mx2),...
  119. disp(Mx3),disp(points),diary off,echo on
  120.  
  121.